Realistic numerical simulations of solar convection and oscillations 

in magnetic regions 



L. Jacoutot^, A. G. Kosovichev^, A. Wray^, & N. N. Mansour^ 

^ Center for Turbulence Research, Stanford University, Stanford, CA 94305 
'^Hansen Experimental Physics Laboratory, Stanford University, Stanford, CA 94305 
^NASA Ames Research Center, Moffett field, CA 94035 



ABSTRACT 

The goal of this research is to investigate how magnetic field affects the dy- 
namics of granular convection and excitation of solar oscillations by means of 
realistic numerical simulations. We have used a 3D, compressible, non-linear 
radiative magnetohydrodynamics code developed at the NASA Ames Research 
Center. This code takes into account several physical phenomena: compressible 
fluid flow in a highly stratifled medium, sub-grid scale turbulence models, radia- 
tive energy transfer between the fluid elements, and a real-gas equation of state. 
We have studied the influence of the magnetic fleld of various strength on the 
convective cells and on the excitation mechanisms of the acoustic oscillations by 
calculating spectral properties of the convective motions and oscillations. The 
results reveal substantial changes of the granulation structure with increased 
magnetic fleld, and a frequency-dependent reduction in the oscillation power in 
a good agreement with solar observations. These simulations suggest that the 
enhanced high-frequency acoustic emission at the boundaries of active region 
("acoustic halo" phenomenon) is caused by the changes of the spatial-temporal 
spectrum of the turbulent convection in magnetic fleld, resulting in turbulent 
motions of smaller scales and higher frequencies than in quiet Sun regions. 

Subject headings: convection — methods: numerical — Sun: oscillations 



Introduction 



Observations of solar oscillations h ave revealed that their pro perties change significantly 



in magnetic regions (iBrown et al.lll992l ). Using SOHO/MDI data lHindman fc Brown! (119981 ) 
found that the power of Doppler- velocity oscillations with frequencies less than 5.2 mHz 
(which corresponds to the acoustic cut-off frequency) decreases with field strength, but the 
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oscillations at higher frequencies become stronger. The regions of enhanced high-frequency 
oscillation power are usually observed at the boundaries of active regions and are sometime s 
referred as "halos" . Similar effect was found for local oscillation modes by lHowe et al.l (120041 ). 
who also found that in magnetic regions the line width of acoustic modes increases at low 
frequencies and decreases at high frequencies. The interpretation of these results is not clear. 
These variations can be caused by the interaction of acoustic waves with magnetic field and 
also by changes in the properties of solar convection and excitation mechanism of acoustic 
waves. It is well-known that magnetic field inhibits convection, thus causing the reduction 
in the power of acoustic waves produced by the turbulent convective motions. However, the 
details of this process are unknown. We address this problem by using realistic numerical 
simulations of solar con vection in the presence o f magnetic field. The realistic numerical 
simulations pioneered by lStein &: Nordlundl (120011 ) have provided an important insight into 
the excitation mechanism of oscillations on the Sun and solar-type stars, showing that the 
acoustic modes by the work of tu rbulent pressure (Reynolds stresses) and nonadiabatic 
gas pressure (entropy) fluctuations (IStein et a/.ll2004r). Of course, the numerical simulation 
must rely on sub-grid scale models of turbulence. IJacoutot et al.l (120081) investigated various 
turbulence models and showed that the dynamic model of iMoin et al\ (jl99ll ) provides the 
best agreement with solar observations. In this paper, we include magnetic field in the 
realistic simulations and investigate changes in the physical properties of convection and 
oscillations. 



2. Numerical model 

We use a 3-D, compressible, non-linear radiative- magnetohydrodynamics code developed 
by Dr. A. Wray for simulating the upper solar convection zone and lower atmosphere. This 
code takes into account several physical phenomena: compressible fluid flow in a highly 
stratified medium, radiative energy transfer between the fluid elements, a real-gas equation 
of state, and magnetic effects. Direct numerical simulation (DNS) of the high Reynolds- 
number turbulent motions on the Sun is not achievable. It is not possible to resolve all the 
motion scales even with the most advanced computers. The large-eddy simulation (LES) 
method allows overcoming the Reynolds-number limits possible for DNS by modeling the 
effects of the smallest turbulent scales. Thus, LES method is used in the present calculations. 

The equations we solve are the grid-cell averaged (henceforth called "averages" ) conser- 
vations of mass ([1]), momentum ([2]), energy ([3]), and magnetic flux (jl]): 
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where p is the averaged mass density, Ui is the Favre-averaged velocity, Bi is the magnetic 
field, and E is the averaged total energy density E = ^pUiUi+ pe+ p(f)+ -^BiBi, where (p is the 
gravitational potential and e is the Favre-averaged internal energy density per unit mass. 
F™*^ is the radiative flux, which is calculated by solving the radiative transfer equation, 
and Pij is the averaged stress tensor P^j = {p + 2puk^kl'i + ^B^Bk) Sij — p {uij + Uj^i) — 
■^BiBj, where p is the vis cosity. The gas pr essure p is a function of e and p through a 
tabulated equation of state ( Rogers et al. 19961 ): Tij is the Reynolds stress, k is the molecular 



thermal conductivity, kt is the turbulent thermal conductivity, a is the molecular electrical 
conductivity, and is the turbulent electrical conductivity. 

We have carried out simulations using th e most widely used subgrid-scale models (SGS): 
the Smagorinsky model (jSmagorinskyI Il963l ) . The Reynolds stresses Tij are modeled by the 
usual eddy viscosity assumption using the large-scale stress tensor Sij = {uij + Uj^i)/2: 



Tij 



-2vt S, 
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(5) 



The eddy viscosity ut is deflned by the Smagorinsky model (jSmagorinskyI Il963l ). ut = 
CsA'^\S\; a nd the trace of th e subgrid Reynolds stress T^k is modeled using the Yoshizawa's 
expression ( Yoshizawa Ill986 ): Tkk = SCgA^ISp, where 15*1 = sj2SijSij and A = (AxAyAz)^/^, 
Ax, Ay, and Az are the grid step sizes. 

The application of the eddy viscosity, pt, and the subgri d kinetic energy, Tkk, results in 
the compressible Smagorinsky formulation (IMoin et a/.lll99l[ ): 



Ti, = -2CsA''\S\{Si, - Uk^kS^jm + 2CcA2|5|25,,/3, 



(6) 



where Cs is the classical Smagorinsky coefficient as used in incompressible flow modeling; 
and Cc is a coefficient associated with the trace of the subgrid Reynolds stress (which is 
absent in the incompressible flow formulation). The two parameters must be specifled in 
some way. For the original Smagorinsky model, constant values are used. In the present 
work, we have chosen Cs = 0.2 and Cc = 0.1. 



The turbulent Prandtl number was taken as unity to set kt- The molecular viscosity p 
and thermal conductivity n were neglected as their solar values are exceedingly small. 
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The turbulent electrical conductiv ity ar is calculated by using the extension of the 
Smagorinsky model to the MHD case (ITheobald et a/.l 119941 ). with the flow velocity shear 
replaced by a shear of the magnetic fleld (which is the electric current): o"t = C^A^IJI with 
ji = ^iikBk.i is the resolved elect ric current; eijk is the Levi-Civita symbol; Cb was taken 
equal to 1. (ITheobald et a/.l 1 1 9941 ). 

In all the following sections, we simulate the upper layers of the convection zone using 
66x66x66 grid cells. The region extends 6x6 Mm horizontally and from 5.5 Mm below the 
visible surface to 0.5 Mm above the surface. The initial magnetic fleld is imposed as a vertical 
uniform compo nent on a snap s hot o f the preexisting hydrodynamic convection, calculated 
as described in lJacoutot et al.l (120081 ). The initial values of the fleld strength are: 100, 300, 
600, and 1200 Gauss. The weak fleld of 100 G has almost no effect. Thus, we present here 
the three other cases. 



3. Granulation structure 

The flrst step was to investigate the influence of the initial magnetic fleld on the granular 
structure of convection. Figure [1] shows temperature, vertical velocity, and vertical magnetic 
fleld distributions at the visible surface for different magnetic flelds. We observe that the size 
of granules decreases as the initial magnetic fleld increases. Without magnetic fleld the mean 
size of granules is about 2 Mm, and it is less than 0.75 Mm for 1200 Gauss. In addition, 
the temperature in the granules becomes higher as the initial magnetic fleld increases. We 
can also note that the down-flow in the intergranular lanes is weaker for high magnetic 
flelds. Figure [21 which shows the vertical velocity energy as a function of the horizontal 
wavelength, k^, measured in degrees of spherical harmonics, / = khR^, demonstrates this 
feature. The vertical velocity energy is decreased by a factor of 10^ for small values I, which 
correspond to the larger granules. We can also notice that the magnetic fleld is swept into 
the intergranular lanes although the magnetic fleld is seeded uniformly. T his characteristic of 



solar magnetoconvection is well-known and has already been presented by lStein &: Nordlund 



(120021 ). 



4. Kinetic energy and oscillation power of radial modes 

We then studied how the kinetic energy is dissipated for different initial magnetic fleld. 
To do this, we calculated the oscillation power spectra of radial modes. Those modes are 
extracted by horizontally averaging the vertical velocity and Fourier-transforming in time. 
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These results are obtained from simulations of 60 hours of solar time using instantaneous 
snapshots saved every 30 seconds. Five oscillation modes can be clearly seen as sharp peaks 
in the spectra of the horizontally averaged, depth-integrated kinetic energy (Fig. [3]) and 
the power spectrum of the vertical velocity oscillations at the solar surface (Fig. Sj). The 
smallest resonant mode frequency is 2.07 mHz. This mode is excited along all the depth. 
The resonant frequencies supported by the computational box are 2.07, 3.03, 4.10, 5.23, 
6.52 mHz. In addition, several br oad high-frequency peaks at 6-12 mHz corresponding to 
pseudo-modes (IKumar fc Lulll99ll ) can be identified. 



It is clear that the magnetic field significantly affects the kinetic energy spectrum. The 
amplitude of the excited modes decreases as the initial magnetic field increases, as well as the 
total oscillation power (Fig. E]). The power distribution is shifted towards higher frequencies 
with the increase of the field strength. It is particularly interesting that the amplitude of 
the pseudo-modes increases with field strength and reaches maximum at Bzo = 600 G. This 
may exp lain the effect of enhanced high-frequency em i ssion ("acoustic halq s") around active 
regions (IBrown et al.l Il992l : iHindman fc Brownl 1 19981 : iJain fc Haberl |2002| ) . The enhanced 
emission at frequencies 5-7 mHz appears at the boundaries of active regions where magnetic 
field is moderate. This corresponds to the simulations results: the pseudo-mode amplitude 
is high for 600 G field and diminishes at 1200 G. 

The simulations also show the enhanced spectral power of the convective background 
at high frequencies for models with magnetic field, forming plateaux aX u > Q mHz. This 
leads to the idea that the acoustic halos are caused by enhanced high-frequency turbulent 
convective motions in the presence of moderate magnetic field. This is consistent with the 
decreased granular size in magnetic regions, described in the previous section. The smaller 
scale convection naturally has higher frequencies and, thus, generates more higher frequency 
acoustic waves than convection without magnetic field. When the field is very strong the 
sound generation decreases because of suppression of convective motions of all scales. This 
probably explains why the acoustic halos are observed in regions of moderate magnetic field 
strength at the boundaries of active regions. 

In addition, the simulation results show that the modal lines in the oscillation power 
spectrum (Fig.H]) becomes broader for the mode of 2.07 and 3.03 mHz (the first mode almost 
disappears at Bzo = 1200 G), but they are more narrow for the r nodes of 4.10 , 5.23 , 6.52 
mHz. This qualitatively corresponds to the observational result of iHowe et al.l (120041 ). We 
plan to present a quantitative analysis in a future publication. 
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5. Calculation of the p-mode work integrand 



The dominant driving of the p-modes comes from the interaction of the non-adiabatic, 

incoh erent pressure flu ctuations with th e coherent mode displacement (INordlund fc SteinI 

200lh Previous work of lBalmforthI Jl992h , IColdreich eTail Jl994l ) , iKumail Jl994[ ). lNordlund fc Stein 
(119981 ) showed that turbulent motions stochastically excite the resonant modes via Reynolds 
stresses (turbulent pressure) an d entropy fluctuat i ons (gas pressure ) . Many analytical formu- 
l ation s have been proposed by iBalmforthI (119921 ): iGoldreich et al\ (119941 ): ISamadi fc Goupil 
(12001! ) ■ The problem is that these models introduce free parameters related to the choice 
of the turbulent medium model. The advantage of 3-D realistic numerical simulations 
comes from the possibili ty of calculating all the quantities related to turbulent convection. 
Nordlund &: SteinI (120011 ) presented a formalism for analyzing the interaction of convection 
with purely radial oscillations. This formalism seems to be the most accurate because it 
accounts for phase relations betw een pressure fluct uations (both turbulent and gas) and the 
mode compression factor dC,^/dr ( IStein et a/.ll2004l ). In this section the mode excitation rate 
is calculated using the same method, but we add the contribution of the magnetic pressure 
-^BkBk in the calculation of the total pressure. The rate of energy input to the modes per 
unit surface area {erg.cm~'^.s~^) is 



A<E^> _ u;^\ldr6P:{d^Jdr 
At SAuE,., 



(7) 



where 6P* is the Fourier transform of the non-adiabatic total pressure; 6 in front of the 
pressure means that one computes so -called 'pseudo-Lagrangi an' fluctuations relative to a 
flxed mass radial coordinate system (INordlund fc Steinll200ll ): Au is the frequency inter- 
val for the Fourier transform; is the mode displacement for the radial mode of angular 
fr equency uo. It is obtained fr o m the eigenmode calculations for a standard solar model 
of IChristensen-Dalsgaard et al\ (119961 ). His spherically symmetric model S gives 35 radial 
modes that provide much denser frequency spectrum in comparison with the three resonant 
modes obtained within the simulation box. E^^ is the mode energy per unit surface area 
{erg.cm~'^) deflned as: 



-UJ 

2 

The distributions of the integrand of the work integ ral as a function of dep th and 
frequency (Figure E]) are similar to the results obtained by lStein fc Nordlundl (120011 ). Most 
driving is concentrated between the surface and 500 km depth at around 3-4 mHz. We can 
see that the excitation becomes weaker for high initial magnetic fleld. The contribution of 
the different pressures (gas, turbulent and magnetic) is presented in Fig. [71 We see that 
the dominant driving comes from the interaction of the coherent mode displacement with 
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the non-adiabatic, incoherent gas pressure fluctuations. The contribution of the turbulent 
pressure shows similar shape than that of the gas pressure, but it is weaker. The magnetic 
pressure plays a weak role in the excitation of the p-modcs for 300 Gauss. This contribution 
becomes higher with increased magnetic field. The total work integrand with increased 
magnetic field becomes localized closer to the surface and at the mode frequencies. 

6. Conclusion 

We have carried out realistic simulations of solar convection and oscillation in the pres- 
ence of magnetic field of various strength. Initially, the magnetic field was vertical and 
uniformly distributed. The results reproduced several phenomena observed in solar mag- 
netic regions. 

In particular, the results confirm that the spatial scale of granulation substantially 
decreases with the magnetic field strength. Magnetic field is swept in the intergranular 
lanes, and the vertical downdraft motions in these lanes are suppressed. This results is a 
decrease in the excitation power. The oscillation power in the presence of magnetic field is 
shifted towards higher frequencies, also increasing the amplitude of pseudo-modes above the 
acoustic cut-off frequency. At a moderate field strength of ~ 600 G the power of the high- 
frequency oscillations reaches a maximum. This corresponds to the phenomenon of " acoustic 
halo" observed in the range of 5-7 mHz at the boundaries of active regions. The reason of 
this phenomenon is probably in the change of the spatial-temporal spectrum of the turbulent 
convection in magnetic field. The convective motions becomes smaller in spatial scale and 
faster in the presence of magnetic field, and this causes the changes in the oscillations excited 
by these motions. 

The power spectra show a significant increase of the convective power ("plateau") at 
high frequencies. The oscillation modal lines becomes broader and smaller in amplitude 
for the low-frequencies modes, but more narrow at higher frequencies. Qualitatively, this 
corresponds to the observations, but requires a more detailed quantitative study. 

The calculations of the work integrand show that the excitation mechanism in magnetic 
regions remains the same as in the quiet Sun. The oscillations are excited by fiuctuations 
of Reynolds stresses and entropy with magnetic forces playing a minor role. With increased 
magnetic field the work integrand becomes more concentrated in the near-surface layers and 
at the mode frequencies. 

In general, our simulations lead to a conclusion that the observed changes of solar 
oscillations in magnetic regions on the Sun are mostly caused by changes in the spatial- 
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temporal spectrum of convective motions, which are the source of the oscillations. 
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Fig. 1. — Temperature {K), vertical velocity {cm/s) and vertical magnetic field (Gauss) 
distributions at the visible surface for different initial vertical magnetic fields. 
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Fig. 2. — Spatial energy of the vertical velocity as a function of the spherical harmonic degree 
/ at the visible surface for different initial vertical magnetic fields. 
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Fig. 3. — Kinetic energy spectra integrated over depth for different initial vertical magnetic 
fields. 
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Fig. 4. — Oscillation power spectra at the visible surface for different initial vertical magnetic 
fields. 
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Fig. 5. — Total oscillation power at the visible surface as a function of initial magnetic field 
strength. 
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Fig. 6. — Logarithm of the total work integrand (eq. ([7j) in units of erg cm as a 

function of depth and frequency for different initial vertical magnetic field. 



-14- 




Depth (Mm) Depth (Mm) Depth (Mm) 



Fig. 7. — Gas, turbulent and magnetic pressure contributions in the work integrand (in units 
of ergcm^'^ ■5^^), as a function of depth and frequency for an initial vertical magnetic field 
of 300 Gauss. 



